setwd("D:/R/class3/class6")
merged_data <- read.csv("result.csv")
rownames(merged_data) <- merged_data$X
merged_data <- merged_data[,-1]

# 将0值标记为NA
pos = merged_data == 0
merged_data[pos] <- NA

# 对数化
merged_data_log2 <- log2(merged_data)

# 替换列名中的数字为字母
group <- gsub("\\d","",colnames(merged_data))

# 存储p值
row_names <- rownames(merged_data)
p_value <- data.frame(P=rep(NA,length(row_names)), row.names=row_names)

# 逐行进行方差分析
for (i in 1:nrow(merged_data)){
  pro <- merged_data_log2[i,]
  pro <- as.vector(t(pro))
  
  # 根据组的名称提取不同的组数据
  posctl <- group == "ctl"
  ctl <- pro[posctl]
  posad <- group == "ad"
  ad <- pro[posad]
  posasym <- group == "asym"
  asym <- pro[posasym]
  
  # 计算每个组的有效数据数量
  ctlnum = sum(!is.na(ctl))
  adnum = sum(!is.na(ad))
  asymnum = sum(!is.na(asym))
  
  # 检查每个组的样本数量
  if (ctlnum > 1 && adnum > 1 && asymnum > 1) {
    # 如果每个组的样本数量足够，则进行方差分析
    anovaresult <- oneway.test(pro ~ group)
    p_value[i,] <- anovaresult$p.value
  } else {
    # 如果某个组的数据不足，则将p值设为NA
    p_value[i,] <- NA
  }
}

# 查看结果
head(p_value)





#题目二
# 安装并加载clusterProfiler包
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
library(clusterProfiler)
library(org.Hs.eg.db)
prostat1<-prostat[prostat$P < 0.05, ]
entrez_gene_ids <- prostat1$ID  # 根据实际的基因ID替换

# 进行GO富集分析
go_enrich <- enrichGO(
  gene = entrez_gene_ids,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",  # 使用符号进行注释
  ont = "ALL",  # 可选择 "BP", "MF", "CC", 或 "ALL"
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

# 查看富集分析结果
head(go_enrich)

# 绘制气泡图
library(enrichplot)
dotplot(go_enrich)
